clc;

% == 1. 导入数据 ==
% 读取 CSV 文件，假设第一列为时间，第二列为信号值
fileName = '30kA0.csv';  % 替换为你的文件名

% 检查是否有标题行（使用 'NumHeaderLines' 来跳过首行标题）
% disp('读取数据...');
data = readmatrix(fileName, 'NumHeaderLines', 1);  % 自动跳过第一行文本

% == 2. 提取时间和信号列 ==
time = data(:, 1);       % 时间列（单位秒）
signal = data(:, 4);     % 原始信号列

% == 3. 确定采样频率 Fs 和归一化截止频率 Fn ==
dt = mean(diff(time));  % 时间间隔均值（假设采样均匀）
Fs = 1 / dt;            % 采样频率 (Hz)

% 奈奎斯特频率
Fn = Fs / 2;

% == 4. 滤波器设计：5阶低通 Butterworth 滤波器，截止频率设为有效频率的 0.1 倍 ==
order = 5;              % 滤波器阶数
Fc = 0.18 * Fn;          % 截止频率（例子中为有效频率的10%）

% 归一化截止频率
Wn = Fc / Fn;

% 设计Butterworth数字滤波器 数字低通 Butterworth 滤波器系数
[b, a] = butter(order, Wn);  % 直接设计数字滤波器

% 或者换成 FIR 滤波器（更平坦滤波特性，如移动平均类）:
% N = 30;
% cutoffFreq = 0.02;  % 归一化截至频率 [0, 1] 对应到 [0, Fs/2]
% lowPassFIR = fir1(N, cutoffFreq);
% filteredSignal = filter(lowPassFIR, 1, signal);

% 将滤波后的信号相位对齐（常用方法避免信号延迟），使用 filtfilt
filteredSignal = filtfilt(b, a, signal);  % 零相位滤波，支持在非实时数据中前后过滤

% == 5. 画图显示原始与滤波后信号 ==
figure;

% subplot(2,1,1);
% plot(time, signal, 'b');
% title('原始信号');
% xlabel('时间 [s]');
% ylabel('幅值');
% grid on;

% subplot(2,1,2);
% plot(time, filteredSignal, 'r');
% title('低通 Butterworth 滤波后的信号');
% xlabel('时间 [s]');
% ylabel('幅值');
% grid on;

% 或者以下方式：将原图和滤波放在同一个坐标轴对比
plot(time, signal, 'b', time, filteredSignal, 'r--');
legend('原始信号', '滤波后信号');
title('滤波前后波形对比');
xlabel('时间 [s]');
ylabel('幅值');
grid on;

% 使用Savitzky-Golay平滑+导数滤波（Summit谱幂滑波序极惯用）  
% ppmRange = Savitzky 修正点 (SF)
filteredSignal = sgolayfilt(filteredSignal, 3, 9);  % 3阶，拟合导数
% slopeNow = slope(zeroCrossingIndices) % 得到每个过零点对应梯度

% == 假设已经导入并滤波后的数据：time 和 filteredSignal
% 示例中 calculatedFilteredSignal = filteredSignal

% 1. 检测零交叉区域的索引
signChanges = diff(sign(filteredSignal(1:1000000,1))) ~= 0;             % 判断信号符号是否改变
zeroCrossingIndices = find(signChanges);                    % 得到前一时刻并找到每次符号变点

time_new = time(1:1000000,1);
% 2. 对每个过零区间进行线性插值处理，得到过零时刻，同时估计斜率
numCrossings = length(zeroCrossingIndices);
zeroCrossingTime = zeros(numCrossings, 1);           % 存储过零点的时间
zeroCrossingSlope = zeros(numCrossings, 1);          % 存储对应的斜率

for i = 1:numCrossings
    idx = zeroCrossingIndices(i);  % 过零起始下标
    t1 = time_new(idx);
    t2 = time_new(idx+1);
    y1 = filteredSignal(idx);
    y2 = filteredSignal(idx+1);
    dt = t2 - t1;
    if dt == 0
        continue;   % 跳过分母为零的点
    end
    % 线性插值计算过零时刻 tz
    if y2 ~= y1
        tz = t1 - y1 * (t2 - t1) / (y2 - y1);
        zeroCrossingTime(i) = tz;
        
        % 计算过零点附近的局部斜率（数值梯度法）
        ts = time_new(idx:idx+1);
        ys = filteredSignal(idx:idx+1);
        slope = (ys(2) - ys(1)) / (ts(2) - ts(1));
        zeroCrossingSlope(i) = slope;
    end
end

% 去除未被赋值的0
zeroCrossingTime = zeroCrossingTime(zeroCrossingTime ~= 0);
zeroCrossingSlope = zeroCrossingSlope(zeroCrossingSlope ~= 0)/1e6;

% 3. 打印&可视化
disp('过零点斜率：');
disp(zeroCrossingSlope);

% 可视化：画出波形和过零点
figure;
plot(time_new, filteredSignal(1:1000000,1), 'b', 'DisplayName', '滤波后信号');
hold on;

% 用红点标出过零点
plot(zeroCrossingTime, zeros(size(zeroCrossingTime)), 'ro', 'MarkerFaceColor', 'r');
legend('滤波后信号', '估计过零点','FontName','黑体','FontSize',10);
% 设置坐标轴范围
axis([-0.002 0.01 -10000 30000]);  % 
% 标注斜率（可选）
for i = 1:length(zeroCrossingTime)
    % 简化标注策略：仅在屏幕上稍作显示
    text(zeroCrossingTime(i), 0, sprintf(' %0.2f', zeroCrossingSlope(i)), ...
         'VerticalAlignment','bottom', ...
         'HorizontalAlignment','right',...
         'FontSize', 10,'Color','red');
end

xlabel('time (s)');
ylabel('Current /A');

set(gca,'FontName','Times New Roman','FontSize',18);
title('过零点及斜率检测','FontName','黑体','FontSize',18);
grid on;
hold off;

